Universality of light thermalization in multimoded nonlinear optical systems

Recent experimental studies in heavily multimoded nonlinear optical systems have demonstrated that the optical power evolves towards a Rayleigh–Jeans (RJ) equilibrium state. To interpret these results, the notion of wave turbulence founded on four-wave mixing models has been invoked. Quite recently, a different paradigm for dealing with this class of problems has emerged based on thermodynamic principles. In this formalism, the RJ distribution arises solely because of ergodicity. This suggests that the RJ distribution has a more general origin than was earlier thought. Here, we verify this universality hypothesis by investigating various nonlinear light-matter coupling effects in physically accessible multimode platforms. In all cases, we find that the system evolves towards a RJ equilibrium—even when the wave-mixing paradigm completely fails. These observations, not only support a thermodynamic/probabilistic interpretation of these results, but also provide the foundations to expand this thermodynamic formalism along other major disciplines in physics.

Generally, for an arbitrary nonlinear photonic lattice, the temperature T and chemical potential µ after reaching the equilibrium can be numerically predicted by the initial condition P, U , M and the spectrum ε j [1, 2]. A recent work [3] shows that, under large M condition, explicit formulas of T and µ can be obtained for a one-dimensional photonic array with uniform nearest-neighbour coupling κ, and that is (1.1) (1.2) Rewriting Eqs. (1.1) and (1.2) can get (1.4) Table 1 list the values of P and U used in the simulations in Fig. 2  The governing equation for cascade χ (2) process is where a m and b m is the field amplitudes associated with the fundamental frequency and its second-harmonic on the mth site, respectively, and ∆ is the phase mismatch. The associated Hamiltonian can be By introducing the generalized coordinates and momenta coordinates (a m , π m , b m , τ m ), where It is easy to check that this Hamiltonian can recover the dynamic equation described by Eq.
(2.1) by following (2.4) The Hamiltonian H does not explicitly depend on z ( ∂H ∂z = 0), and therefore H is conserved. Here the Poisson bracket is defined as: where q i and p i are the generalized coordinates and canonical momenta.
The second conservation law in this system is . This can be checked by the Poisson bracket (2.7) Therefore

Supplementary Note 3: Conservation laws in nonlinear optomechanical cavities
The governing equation for a nonlinear optomechanical cavity array is where a m and b m stands for the optical field and the mechanical oscillation amplitude in cavity m, respectively. Here the parameter Ω represents the frequency detuning. The Hamiltonian is The coordinates and canonical momenta are (a m , π m , b m , τ m ) and π m = −ia * m , τ m = −ib * m . Then the Hamiltonian can be written as On the other hand, P b = M m=1 |b m | 2 = M m=1 ib m τ m for the mechanical oscillators is not conserved, since [P b , H] ̸ = 0 as follows: (3.5)

Supplementary Note 4: Smooth but nonanalytic function
A smooth non-analytic function f 1 (x) can be generated numerically by using the method described in Ref. [4]. One example of such type of functions is shown in Fig. S1a. This function is nowhere analytic in its domain since its Taylor series expansion does not converge. This is illustrated in Fig. S1b, where we observe that the magnitude of the first, second and third derivatives of f 1 (x) grows from 50 to the order of 10 6 . In our study of light thermalization, we considered only a part of this f 1 (see red line in Fig. S1a) and performed the following transformation to obtain F 1 (x) shown in Fig. S1c): FIG. S1. a The smooth everywhere but nowhere analytic function f 1 (x) constructed by Fourier series with random coefficients of root-exponentially decreasing amplitudes. b The first, second and third derivatives of function f 1 (x). c In our simulation, we use a function F 1 (x), which is a part of f 1 (x) after a transformation according to Eq. (4.1).

Supplementary Note 5: Field intensity distribution in the saturable nonlinearity
In the simulation of a saturable nonlinearity in Fig. 3f of the main text, to ensure that the Laurent series expansion is valid (i.e., |a m | 2 > 1), the values of the local intensities |a m | 2 have been monitored, and their distribution is illustrated in Fig. S2. As we can see, only 5% of |a m | 2 locate in the range of [0, 1] and the rest are greater than unity .
FIG. S2. The distribution of the values of |a m | 2 during simulation.

Supplementary Note 6: Thermalization in a SSH lattice
The light thermalization is a general phenomenon, which can happen in a lattice with uniform coupling as discussed in the main text, but also can occur in topological lattice with nonuniform coupling coefficients. In this section, we investigate the optical thermalization in a topological Su-Schrieffer-Heeger(SSH) lattice (Fig. S3a). The coupling coefficients are taken κ 1 = 4/3, κ 2 = 2/3, and the coresponding spectrum ε j is plotted in Fig. S3b. In this lattice, an input excitation |c j | 2 = 0.05(ε j + 2) (dashed line in Fig. S3c) leads to P = −10.9, which in turn predicts the RJ distributation with T = 0.13 and µ = −2.4, which is consitent with the numerical simulation result as shown in Fig. S3c.
Moreover, a non-Hermitian SSH model with asymmetric couplings between nearest neighbors was reported in [5] recently.